Probability Density Function of Longitudinal Velocity Increment 
in Homogeneous Turbulence 
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Two conditional averages for the longitudinal velocity increment Ur in the simulated turbulence 
are estimated; h(ur) is the average of the difference of the Laplacian of the velocity field with the 
Ur value fixed, while g{ur) is the corresponding one of the square of the difference of the velocity 
gradient. Based on the physical argument, we suggest the fitting formulae for h and g, which quite 
satisfactorily reproduced to the 512^ DNS data. The predicted PDF is characterized as (1) the 
Gaussian distribution for smaller amplitudes, (2) the exponential distribution for larger ones, and 
(3) a prefactor before the exponential function for the intermediate ones. 

PACS numbers: 47.27.Gs 



0^ 



> 

o 

o 
m 
o 

00 

ON . 

"a : 
I 

o 



X 
S3 



Probability density function (PDF) is one of the most 
interesting subjects in turbulence The intermittent 
effect often referred to is reflected in the deviation of 
the PDF from the Gaussianity. In particular the PDF 
for longitudinal velocity increments has been the main 
concern in turbulence |2|-||]. For large separation r the 
PDF is Gaussian and it becomes deviated more as r de- 
creases. The PDF of its derivatives is even flatter than 
the exponential distribution In the present letter 

we propose a simple form of PDF, which is consistent 
with the observation so far. 

The conditional average method Jll|-p^ is very useful 
to understand the nature of turbulence. In a previous 



paper |14| the method is applied to the study of the pas- 
sive scalar field advected by Navier-Stokes turbulence. In 
the present letter the method is applied to the analysis 
of statistics of the longitudinal velocity difference at two 
points separated by r: ?7r(x) = ?7(x-|-rx) — C/(x), where 
t/(x) is the X component of the velocity field. The first 
conditional average we study is 



H{Ur) = {VlUr{^)\Ur), 



(1) 



where the average is taken with a value of Ur fixed. The 
second average is 



G(C/.) = (|VxC/r(x)nc/,). 



(2) 



For later convenience U , H, and G are made dimen- 
sionless as 



9{Ur) = -^G{Ur) 

Mr) 



Mr) 



(3) 



where v is the kinematic viscosity, S2{r) = {Ur}, and 
J2(r) = -2z/(|VxC/.(x)|2). 



A useful relation concerning the PDF P{u) of u is de- 
rived for a statistically homogeneous system in terms of 
h and g as follows [nUTlJTll: 



P{u)g{u) = Cexp 



" h{x) 



9{x) 



dx 



(4) 



where C is a numerical constant. This relation was nu- 
merically confirmed using our DNS data. 

To begin with, let us derive possible forms of h(u) and 
g{ii). The equation for Ur is formally written as 



^+T{Ur)^vV^Ur 

at 



(5) 



where T{Ur) is the inertial term containing the pressure 
effect. Multiplying the above equation by Ur we take the 
conditional average with Ur fixed: 



l/dUl 
2 \ dt 



Ur 



Ur{T{Ur)\Ur) = vUr (V^ Ur\Ur) . (6) 



In the quasi-steady state the first term on the left hand 
side can be neglected. Then 



vUrH{Ur) = Ur{T{Ur)\Ur), 



(7) 



which implies that vUrH{Ur) balances the energy trans- 
fer rate in the inertial region. It would take the form 



Ur{T{Ur)\Ur) ^ ^, 



(8) 



where the relaxation time would be given by the eddy 
turnover time r/Ur, then Ur{T{Ur)\Ur) ~ Ur/r. Hence 



r 



(9) 
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In the inertial range, the eddy turnover time r,. is given 
by the Kolmogorov scahng form ~ r^/'^e"^/'^ (e is the 
average dissipation rate), independent of Ur, which im- 
phes that 



vH{Ur) ^ e 3 r ^Ur ccUr 



(10) 



The scahng (|^) wih be dominant for large values of Ur 
and the Kolmogorov scaling (10) is for small values of 
Ur- Since mechanisms coexist, a combined form 



h{u) = cqu + c\v^ 



(11) 



is expected in the nondimcnsional form. We have at- 
tempted a linear fitting for small values of u and a 
quadratic one for large values of u, but the fitting ( [ll| ) 
over the entire region is more satisfactory. As the 
Reynolds number increases, the range proportional to iP' 
becomes wider. 

Next we consider G{Ur)- Multiplying (||) by Ur and 
taking a certain average, we have 

{UrT{Ur)) ^ ■ (UrVUr)) ~ iy{\yUr\^) . (12) 

The first term on the right hand side vanishes because of 
the statistical homogeneity. Then 



Ur{T{Ur)\Ur)P{Ur)dUr = -V / G{Ur)P{Ur)dUr 



Thus G{Ur) is the conditional dissipation rate defined af- 
ter the spatial average being taken. Consider the extreme 
case where the spatial mixing is strong. Then the con- 
ditional dissipation rate G{Ur) should be independent of 
Ur- Another limit is that G{Ur) must be proportional to 
the energy density C/^. Hence one may expect a nondi- 
mensionalized form g{u) = ap + 02^^. (Since the spatial 
transfer works to diminish the scaling, 02 is expected 
to decrease as r increases.) Such a fitting is considerably 
good, but not excellent. If one adds a linear term as 



g{u) = flo + aiu + 02^^, 



(13) 



the agreement is remarkable as we see later. A linear 
term is necessary for a smooth transition from the homo- 
geneous to the intermittent dissipation. 

Now we turn to h{u) and g{u) based on the DNS data 
and the discussion of the PDF. The simulations were car- 
ried out with resolutions of 256'^ and 512'^ mesh points 
using the NAL Numerical Wind Tunnel for a couple of 
Reynolds numbers ||l^. Since most of the results ob- 
tained are very similar to each other, we present only 
the results of the run of 512^. The turbulence decays 
from the initial spectrum fc'*exp(— 2fc^) with no forcing 
and the data were processed when the energy dissipa- 
tion rate reaches the largest value. At that time the 
Taylor microscale Reynolds number is 120. The inertial 
region in which (?7^)/(er) is nearly constant is located 
in 14 < r/77 < 51, 77 being the Kolmogorov microscale. 



The three functions, h{u), g{u), and P{u)g{u) are esti- 
mated by using the simulation data which are presented 
in the following. 

(i) h(u). The curves of h{u) against u are depicted for 
various values of r in Fig. 0. The interesting points are 
noticed by an inspection, (a) All the curves for the iner- 
tial range separations fall on an almost single line, which 
indicates that the conditional expectation of the energy 
transfer rate is independent of the scale, (b) h{u) is not 
antisymmetric; in fact the absolute magnitude is larger 
for large values of u in the positive region than in the 
negative region. This tendency is obvious, because the 
longitudinal component is similar to Burgers turbulence 
|]l7| , where the positive slope region [u > 0) works on 
the negative slope to make it steeper, which implies that 
the positive region transfers energy to the negative slope 
side. The excess of h{u) in the positive region repre- 
sents a faster decay of the PDF there, (c) /i(0) is slightly 
negative for the inertial separations; h{Q) = —0.05 for 
r/77 — 20. This is responsible for the skewed PDF. 




FIG. 1. h va u for various values of r. The solid lines in 
the inset represent the fitting for r/rj — 20, the crosses being 
the data points. 



Now we fit the formula (|11|) with the data. Since the 
curves are not symmetric, we make a separate fit for u > 
and u < 0. The fitting is remarkable except for the 
largest values of where the data are scattered because 
of the insufficient number of data. The inset in Fig. |l| is 
an example of the curve for r/rj — 20 fitted in the interval 
\u\ < 8, where the solid line is the fitting function. The 
fitting parameters are as follows: For r/77 — ^ cq — 0.67, 
ci — 0.27 in the region u > 0, and cq = 0.56, ci = —0.20 
in the region u < 0. On the other hand, for r/77 = 20 i^i 
the inertial region, cq = 0.47, ci = 0.33 for u > 0, and 
Co = 0.44, ci = -0.20 for u < 0. 

(ii) g{u). The curves of g{u) against u are given for var- 
ious values of r in Fig. |^. (a) The curves depend on 
r, which is responsible for the scale dependence of the 
PDF. (b) As r increases, g{u) becomes flatter and g{0) 
increases, (c) g{u) fluctuates considerably for large am- 
plitudes; in particular, more conspicuous for large values 
of r. The reliable range is wider for ?i < than for u > 0. 
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The formula (|l3|) is fitted to the data separately for 
the regions (I) u> bi, (II) w < 62, and (III) 61 > m > 62; 
where 61 = 0.3 and &2 = —0.2 for r/77 — 20, for instance. 
(If only two regions m > and m < were chosen, there 
would be a small discrepancy near u — 0.) The fitting 
is quite satisfactory. The solid line in the inset of Fig. |^ 
represents the fitting functions for r/rj = 20, where the 
fitting is processed for —5 < u < —0.2 and 5 > u > 0.3. 

In order to ascertain that the present result is consis- 
tent with the observation [^|J^], which reports the scale 
dependence of the factor Ci/a2, we listed parameters ag, 
ai and a2 for various values of r in Table ^ Since the 
region (II) is wider than the region (I), only parameters 
fitted for the region (II) are given; ag is found to increase 
steadily with r. The coefficient 02 decreases with r in the 
viscous and lower inertial ranges, but the decrease can- 
not be seen for larger values. Note that the determina- 
tion of a2 is a little subtle, because the region exhibiting a 
quadratic dependence of g{u) cannot be appreciated with 
increasing r in the simulations of low Reynolds numbers. 
This scale dependence of 02, however, does not contra- 
dict with the observation According to Ref. ||] 
Ci/a2 ~ r'^ i'^j meaning that 02 ^ ^--0.17 i3g(2aygg jg 
dependent of r in the inertial range; even when the scale 
is doubled, 02 is reduced only by a factor 0.88. Hence 
the variation of 02 may be too small to be observed in 
the limited interval of inertial range. In this situation 
we conclude that 02 decreases as r increases from the 
dissipative to inertial scale. 
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FIG. 2. g vs u for various values of r. The solid lines in 
the inset represent the fitting functions in the range (I) and 
(II) for r/rj — 20, the crosses being the data points. The 
discrepancy for small amplitudes is seen. Those components 
are fitted by another quadratic expression. 

(iii) P{u)g{u). The expression (^ suggests that the com- 
bined form P{u)g{u) is simpler than P{u) itself. Figure^ 
is a collection of P{u)g{u) with varying r, where the nor- 
malization J g{u)P{u)du = I is employed. The Gaussian 
part at small amplitudes and the exponential-like part at 
larger amplitudes are perceptible, although the slope of 
the tail becomes steeper as r increases. The PDF is ob- 
tained by dividing P(u)g{u) by g{u). The inset of Fig. ^ 
is the PDF for r/rj — 5, which exhibits the stretched 



exponential exp[— a;''] with p < 1. Hence the difference 
between the exponential distribution in P{u)g{u) and the 
stretched exponential one in P(u) is due to the prefactor 
l/g{u)\ it is checked directly by using the fitting for- 
mula for r /t] — 5, although the result is not shown here. 
It should be emphasized that the prefactor is important 
only for small and intermediate amplitudes, and its effect 
tends to vanish as u becomes large. 



TABLE I. The fitting parameters ao,— oi and a2 for 



various 


values of r in 


the region 


-5 < 


< -0.2. 
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40 
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ao 


0.40 


0.52 


0.67 


0.75 


0.85 


0.88 
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0.59 


0.60 


0.43 


0.33 


0.18 


0.14 


02 


0.26 


0.10 


0.048 


0.039 


0.054 


0.050 
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FIG. 3. P(u)g(u) VS u for various values of r. The dotted 
line is the normalized Gaussian with g{u) = 1. The inset is 
P(u) for r/?7 = 5, where the stretched exponential can be 
sensed for intermediate amplitudes. 

(4) qiu) = h{u)l g(u). The ratio function q{u) determines 
the skeleton of the PDF. Since 

HiUr)dU, ^ Ur{T{Ur)\Ur) dUr 

G{Ur) vG{Ur) Ur' ^ ' 

it is reasonable to think that q{u) is a ratio of the energy 
transfer with fixed u,. to the spatial averaged dissipation 
rate with fixed Ur- Figure ^ is the curve q{u) for var- 
ious values of r/rj. For r/r] = 5, it looks the shifted 
tanh(u). As r increases the curves are gradually devi- 
ated from tanh(u), which is consistent with the fact that 
g{u) is dominated more weakly by the term a2M^ since 
02 decreases with r. It is imagined, however, that as the 
system size goes to infinity, q{u) approaches a constant 
value for large u, implying the exponential tail for P{u). 

The PDF of the longitudinal velocity increment is 
skewed, and its effect is reflected in the property that q(u) 
is not antisymmetric about the origin, as seen from Fig.^. 
It is interesting to notice, however, that the curves q{u) 
vs u are antisymmetric about a point marked by an ar- 
row. In fact, the inverted curve (a dotted line) about the 
marked point in the inset of Fig. coincides almost with 
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the original one. When the point is denoted as {n,m), 
n ^ m ^ 0.2 for any separation. Under this situation 
the function q{u) can be written as q{u) = f{u — n) + m, 
where f{u) is an antisymmetric function. 
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FIG. 4. q vs u for various values of r. In order to see 
the antisymmetry of q{u) about a point marked by an arrow, 
we depicted the curve for r/i] = 20 (a solid line) and the 
inverted curve (a dotted line) about the point in the inset. 
The difference is not appreciated. 



For vi{r) < Ur < V2{r) it is of exponential type. When r 
increases with the Reynolds number fixed, vi{r) becomes 
of the same order as V2{r) at r — vq. Then for r > vq the 
PDF never becomes exponential. If the Reynolds num- 
ber is small enough, the situation V2{r) < vi{r) occurs 
for any r, implying that the exponential distribution can- 
not be observed, in agreement with the observation that 
the exponential distribution was identified only in high 
Reynolds numbers |^,^. 

When we consider the PDF in a fully developed turbu- 
lence, the following questions arise: what is the PDF of 
smaller scales when the Reynolds number increases? Is it 
the stretched exponential as e~'^'' with p going to zero? 
The present study answers these questions on the basis of 
the physical argument supported by the simulation. The 
answer is that the eventual PDF is the exponential for 
large amplitudes. Only the slope increases with r. The 
stretched exponential is spurious due to the prefactor and 
it vanishes for large amplitudes. 



A conjecture for such an expression is as follows. Let 
us start with an antisymmetric function f{u): 



P{u)g{u) = Cexp 



f{x)dx 



Since the third order moments of u (skewness) must be 
negative because of the energy transfer, we must have 
a skew P{u). {g{u) is approximately symmetric around 
the region of u which we are concerned.) The simplest 
way to do so is to multiply the above expression by e~™". 
However the obtained PDF yields the negative first order 
moment. In order to make the first order moment zero, 
we must shift u by a certain amount n as 



P{u)g{u) = Ce- 



' exp 



f{x)dx 



which agrees with the observation. Unfortunately we 
have not known a reason for the conjecture. 

The present result is summarized for the PDF except 
for the prefactor effect signified by 1/ g{ur)- For simplic- 
ity we ignore its asymmetric effect. First we consider 
the limit of the infinite Reynolds number. In the region 
with small values of Ur the PDF is Gaussian; P(ur) ^ 
Q-a{r)u^^ where a{r) — co(r)/(2ao(r)) is a decreasing 
function of r. For large values of Ur P{ur) — e^''^''^"'', 
where (3(r) — Ci{r)/a2{r) is an increasing function of r. 
The transition value vi{r) is of order of (3{r) / a{r), which 
is an increasing function of r. Around ^ viir) the 
PDF transits from Gaussian to exponential. 

Now we turn to the case of finite Reynolds numbers. I n 
this case large amplitudes of Ur with a fixed r cannot exist 
beyond the cutoff amplitude f 2 {r) , which increases with 
the Reynolds number. Then we have the following PDF. 
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